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The change in the decay rate of an excited atom that is brought about by extinction and thermal- 
radiation effects in a nearby dielectric medium is determined from a quantummechanical model. The 
medium is a collection of randomly distributed thermally-excited spherical scatterers with absorptive 
properties. The modification of the decay rate is described by a set of correction functions for which 
analytical expressions are obtained as sums over contributions from the multipole moments of the 
scatterers. The results for the modified decay rate as a function of the distance between the excited 
atom and the dielectric medium show the influence of absorption, scattering and thermal-radiation 
processes. Some of these processes are found to be mutually counteractive. The changes in the 
decay rate are compared to those following from an effective-medium theory in which the discrete 
scatterers are replaced by a continuum. 

PACS numbers: 42.50.Nn, 42.50.-p, 41.20.Jb 


I. INTRODUCTION 

The emission process of a photon by an excited atom is 
influenced by its environment m- In particular, when 
the excited atom is located in the vicinity of a dielectric 
body, its decay rate will be different from that in vac¬ 
uum. It will depend crucially on the precise properties 
of the dielectric medium for several reasons. Interference 
effects due to the reflection of emitted photons from the 
surface of a dielectric medium will alter the atomic de¬ 
cay rate. Furthermore, if the medium is dispersive and 
lossy, it may absorb photons. This absorption process 
will lead to an enhanced atomic decay. Scattering effects 
due to inhomogeneities in the dielectric medium will also 
change the atomic decay. Finally, if the dielectric has got 
a finite temperature, it will emit thermal radiation which 
will stimulate the emission process. All these effects will 
depend on the shape of the dielectric body and on its 
distance from the atom. 

Detailed studies of the changes in the atomic decay 
rate due to the presence of homogeneous lossy dielectric 
media at zero temperature have been carried out for var¬ 
ious geometries. The simplest case is that of a halfspace 
that is uniformly filled with a dispersive dielectric [IHZl- 
A more complicated conhguration is that of a dielectric 
sphere, which has been studied extensively [8Ul5]. Other 
geometries like cylinders or spheroids have been consid¬ 
ered as well m- 

In a realistic medium extinction of radiation is pro¬ 
duced not only by absorption but also by scattering from 
inhomogeneities. Both these effects will alter the decay 
rate of a nearby atom. When the inhomogeneities are 
densely distributed, multiple-scattering effects greatly 
complicate the analysis. For such systems numerical sim¬ 
ulations have been used to gain insight in modified emis- 
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sion processes [IIHID]. For a dilute distribution of scat¬ 
terers, however, multiple scattering plays a minor role, so 
that analytical methods may be used. In a recent paper 
m we have studied a model in which an excited atom 
is situated near a halfspace that is filled with a dilute set 
of spherical scatterers consisting of absorptive dielectric 
material. For this model detailed results for the change 
in the atomic decay rate as a function of the distance 
between atom and halfspace have been obtained. 

The changes in the atomic decay rate due to thermal 
radiation from a nearby dielectric medium are of a dif¬ 
ferent nature from those discussed above, as they are 
a consequence of stimulated-emission effects. Although 
their origin is different, they are likewise sensitive to the 
geometric details of the dielectric medium. For a uniform 
medium with a flat interface thermal-radiation effects in 
the atomic emission process have been studied in [22ll23) . 

In the present paper we shall treat within a single 
model decay-rate changes due to both extinction effects 
(from absorption and scattering) and to thermal radia¬ 
tion. We shall consider a randomly distributed dilute set 
of spherical scatterers consisting of a lossy dielectric ma¬ 
terial. By choosing the scatterers to be spherical as in [2T] 
we are able to use the framework provided by Mie the¬ 
ory. Taking the set of scatterers to be dilute permits us 
to neglect multiple-scattering effects, as explained above. 
The scatterers will be held at a finite temperature so that 
they will emit thermal radiation. We will choose the set 
of scatterers to be bounded, so that the optical depth of 
the aggregate stays finite. For analytical purposes, the 
shape of the aggregate will be taken to be spherical. The 
model will allow us to treat the effects of reflective inter¬ 
ference, absorption, scattering and thermal radiation on 
an equal footing. 

The paper is organized as follows. In Section [IT] we in¬ 
troduce our model that is based on the damped-polariton 
formulation of absorptive dielectrics. We formulate the 
atomic decay rate in terms of the second moment of the 
electric field and the associated Green function. In Sec- 
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tion III we use addition theorems for vectorial spheri¬ 
cal wave functions to evaluate the correction functions 
that govern the changes in the atomic decay rate in the 
presence of a set of absorptive scatterers at finite tem¬ 
perature. Subsequently, an alternative approach to the 
correction functions is given in Section [TV] It makes use 
of integral representations, which are more suitable for 
numerical evaluation. Some typical examples of the en¬ 
suing graphs for the correction functions are presented 
in Section |V] Furthermore, a comparison is made with 
an effective-medium model in which a uniform dielectric 
replaces the set of discrete scatterers. In a final Sec¬ 
tion |Vl] our results are summarized and discussed. A few 
technical details of our treatment are given in a set of 
Appendices. 


II. ATOMIC EMISSION AND ABSORPTION 
NEAR A DISPERSIVE DIELECTRIC MEDIUM 


A dispersive and absorptive linear dielectric medium 
may be described by a damped-polariton model that has 
been introduced some time ago [21] . In this model damp¬ 
ing effects are taken into account by coupling the polar¬ 
ization density to a bath of harmonic oscillators with a 
continuous range of frequencies. For a uniform dielectric 
medium Fourier expansions can be used to diagonalize 
the model Hamiltonian and to determine its eigenmodes. 
For arbitrarily inhomogeneous media diagonalization can 
likewise be carried out by means of a Green function tech¬ 
nique |2S]. In the latter case diagonalization leads to a 
Hamiltonian of the form 

H = J J duj huj C\r,uj) ■ C{r,uj), (1) 

with annihilation operators C(r,a;) and associated cre¬ 
ation operators depending on position and frequency ar¬ 
guments. These fulfill the standard commutation rela¬ 
tions [C(r,a;), Cl'(r',w')] = I (5(r — r')(5(a; — w'), with I 
the unit tensor. The electric field can be expressed in 
terms of these operators as E(r) = duj E(r, w) -I- h.c., 
with 


E(r,w) = J dr' fE(r, r', w) • C(r', w). 


( 2 ) 


The tensorial coefficient f e reads [25] : 

1/2 

G(r,r',a;-|-zO). 

(3) 

Here e is the complex local (relative) dielectric constant, 
which follows from the parameters of the model. Further¬ 
more, G is the tensorial Green function, which satisfies 
the differential equation 


^ hlmeir'+ iO) 

f£;(r,r ,a;) =--- 

TT£q 


— V X [V X G(r,r',w-|- jO)] 

-I-— e(r, oj + iO) G(r, r', uj + iO) = I S(r — r'). (4) 


The rate of photon emission by an excited atom in the 
vicinity of a dispersive and absorptive dielectric follows 
from the inhomogeneous damped-polariton model in its 
diagonalized form by employing perturbation theory in 
leading order and disregarding transient effects |2S]. In 
the electric-dipole approximation it can be expressed as 
an integral over the second moment of the electric field: 

Ptt 

^ ■ J dJ (E(ra,w)E^(ra,a;')) ■ {9\^J■\e) , 

(5) 

with atomic dipole matrix elements connecting the ex¬ 
cited state (labelled e) and the ground state {g) of the 
atom, with ui the atomic transition frequency and with 
the atomic position. The second moment of the electric 
field is determined by the density matrix of the damped- 
polariton system. 

After insertion of ([^ one encounters the second mo¬ 
ment of the creation and annihilation operators. This 
second moment will be assumed to be isotropic and di¬ 
agonal in both the position and the frequency variables: 

(C(r,a;) C'l'(r',a;')) = [1-f ^(r, w)] I (5(r - r')(5(a; - w'), 

( 6 ) 

with a real and positive function g{r,uj) that vanishes 
when the damped-polariton system, as given by the 
Hamiltonian 0. is in its ground state. The second mo¬ 
ment of the electric field now gets the form 

fe, A 

(E(r,a;)E'f(r,w')) = i5(w - uj ') -j 

neoc^ 

X J dr'G(r, r', w-I-iO) • G (r, r', w-I- jO) 

x[l-|- 5 (r',a;)] Ime(r',w-I-iO), (7) 

with G the transpose of G. The part of the integral that 
is independent of g(r, uj) can be rewritten with the help 
of the optical theorem for the Green function: 

J dr'G{r,r',uj + iO) ■ G (r, r',a;-I-*0) Ime(r',a;-I-lO) = 

c2 

=-^ ImG(r,r,a;-I-lO). (8) 

To evaluate the second part of the integral in 0 one 
needs information on the function g(r', uj) in the inte¬ 
grand. If g{r',uj) stays finite for large |r|', the distant 
regions in the integral may contribute even when Ime 
tends to 0 there. The importance of the far parts of space 
while evaluating the second moment of the electric field 
has been discussed previously m- In fact, the Green 
function G(r, r',a; -I- iO) for a system with a uniform di¬ 
electric constant e at large distances from the origin will 
be proportional to l/'=/|r'| for large |r'|, so that upon 
integration these distant regions yield a finite contribu¬ 
tion even when Ime gets negligibly small. However, if 
g(r',uj) tends to 0 for large |r'|, as is reasonable when 
no incoming fields are present, the contribution of the 







3 


far spatial domain to the integral is no longer important 
when Ime gets small there. In that case the integration 
in the second part of the integral in 0 can be confined to 
positions in space where Im e differs from zero by a physi¬ 
cally relevant amount, as will be indicated by a subscript 
V at the integral. Upon assuming moreover for simplicity 
that both g(r',uj) and Ime(r', w-I-*0) are independent of 
position within U, we arrive at the following expression 
for the second moment of the electric field: 

(E(r, w)E^(r,a;')) = —8{y} — w')-^ ImG(r,r,a; -I- lO) 

TTEqC^ 

huj'^ 

-|-i5(a; — to') - jgi^) [Ime(w -I- iO)] 

TTSqC^ 

X f dr'G{r,r',uj + iO) ■ G (r, r', w-I-lO). (9) 

Jv 

For a dielectric at a finite inverse temperature /3 one may 
insert — 1). 

The atomic decay rate follows by substituting 0 in 
([^. As a result the decay rate is found as the sum of 
two terms, one representing spontaneous emission in the 
presence of an absorbing dielectric at zero temperature, 
and the other a correction describing stimulated emission 
due to the finite-temperature radiation from the dielec¬ 
tric. These two contributions have different properties 
and should be treated separately, as has been remarked 
before Pl[?fl[^- 

In the presence of thermal radiation from the dielectric, 
atomic transitions from the ground state to the excited 
state may occur as well. The ensuing rate of atomic 
absorption is determined by an expression similar to ([^ : 

27r 

• J du}' (E'f(ra,w)E(ra,w')) ’ (gImIs) • 

( 10 ) 

The second moment of the electric held occurring here is 
given by the analogue of (|^ , with no hrst term and with 
a second term that follows by interchanging G and G*. 


III. DECAY NEAR A SPHERE WITH 
ABSORPTIVE SCATTERERS 

We consider a spherical domain with radius R that 
contains a dilute set of non-overlapping spherical scatter- 
ers with radii a R and complex dielectric constant e, 
while the space between the scatterers is empty vacuum. 
The region U in ([^ is given by the union of the interiors 
of the spherical scatterers. The enveloping spherical do¬ 
main with volume dirR^/S will be indicated by Vq. The 
scatterers are assumed to be randomly distributed with 
a uniform average density. An excited atom is located 
outside the spherical domain, at a distance Va > R + a 
from its center, which is chosen as the origin of a spheri¬ 
cal coordinate system. We wish to determine the change 
in the decay rate that is brought about by the presence 


of the set of spherical scatterers. In this section we will 
first concentrate on the effects of scattering for the case of 
cold scatterers. Subsequently, finite-temperature effects 
will be considered as well. 

For cold scatterers the second moment of the electric 
field is given by the first term of ([^. For a single cold 
scatterer the Green function G is the sum of a vacuum 
contribution Gg and a scattering term Gg, as given in 
Appendix The imaginary part of the vacuum Green 
function with coinciding position arguments is propor¬ 
tional to the unit tensor: ImGo(r, r, w-l-iO) = —wl/(67rc), 
whereas the imaginary part of the scattering term Gg for 
coinciding positions follows from (Al). 

For a dilute set of Ng scatterers, with centers at the 
positions r^, the Green function can be approximated as 


G(r, r',uj + iO) = Go(r,r',a; -I- iO) 

-f ^ Gg(r - ri,r' - ri,a; -b iO), (II) 

i 


as follows by suppressing all terms involving multi- 
scatterer configurations in the Foldy-Lax equations or, 
more directly, in the Neumann series for the integral 
equation that determines the Green function [29l [30] . 
Adopting this form of the Green function implies that 
multiple-scattering effects are henceforth neglected, as is 
reasonable for a dilute set. Since the centers of the scat¬ 
tering spheres are located at random positions inside the 
spherical domain with radius R, the configurational av¬ 
erage of the second term of © may be written as an 
integral over Vg: 


'Vo 


dr"Gg(r — r",r'— r",a;-b lO), (12) 


with rig = 3Ng/{ATrR^) the density of the scatterers. 
Upon inserting the expression (AI) for Gg, employing 
the addition theorem ( |A6[ ) for the vector spherical wave 
functions M[^(r), N|^(r) and carrying out the integral 
over r", we get for the imaginary part of the configura¬ 
tionally averaged Green function with coinciding position 
arguments: 


Im (G(r, r,a; -b iO)) = —-^1 -b 47r—ng l{l + 1) 

OTTC <' 


U) 

' —1 

c 


1=1 


^ ^ (Tylff I 1 


Z' = l m'=-l' l"=0 

xRe 


(-*)'+isr {5,V+,,N|,';^,(r)Nj,%(r) 
+<5?+,+z.M(';^,(r)M(y„,(r)} 

+ (-i)'+ii?r{5?+.+r'M|fL(r)M|,%(r) 


I -I 0 


(13) 


Here Ii{R) = /p^drr^[j;(fcr)]^, with k = w/c, results 
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from the integration over the spherical domain. Its ex¬ 
plicit form is 

hiR) = + 

-1(21 + l){Ryk)jiikR)ji+i{kR ), (14) 

as may be checked by differentiating with respect to R 
and employing the recursion relations for the spherical 
Bessel functions m- The Wigner 3j-symbols in 


imply that the three summation variables 1,1', I" satisfy 
the triangular conditions \l' — I"\ < I < I' + I". Further¬ 
more, Sf equals 1 for even I, and 0 for odd I, while d° is 
defined analogously, with even and odd interchanged. 

Substitution of (13) in the first term of ([^ yields an 
expression for the configurational average of the second 
moment of the electric field. From the spherical sym¬ 
metry it follows that the resulting tensor is diagonal in 
a spherical coordinate system with origin at the center 
of the aggregate of scatterers. Upon using the expres¬ 
sions for the vector spherical wave functions in spherical 
coordinates [35] one gets 


;(E(r,cc)Et(r,cc'))) = <5(cc-cc'), 




X {[l + /Fej|(r)] + [1 +/Fc,_L(r)] (e^ee-b , 

(15) 

with the combined configurational and density-matrix 
average indicated by double brackets. Furthermore, e,., 
eg and are unit vectors in spherical coordinates, and 
/ = I Trains is the filling fraction of the spherical do¬ 
main containing the scatterers. The decay-rate correc¬ 
tion functions Fcj and Fc^± are 


In view of ([^ the average decay rate for an excited 
atom in the presence of a randomly distributed set of 
cold spherical scatterers can be written as 

(kc) = Fo + / [F,,||(r,)ro.|| + Fe,i(r„)ro,j.] , (19) 

with the vacuum decay rate Fg = Fo^n + Fg^j. for Fg y = 
uj^\{e\fi-er\g)\'^/{3TTeohc^) and Fg,_L = uj^{\{e\fi-eg\g)\'^ + 
\{e\fj, ■ e,p|(/)p)/(37regfic3). This expression for the decay 
rate in the presence of an aggregate of cold scatterers 
reduces to a simpler form if only a single scatterer is 
present. The latter form is found by taking the limit R —>■ 
0 and ris —t oo with 47ri?3j^^/3 = go that Ug ci^iijii{R) 
gets equal to As a consequence, only the 

terms with in the correction functions survive 

in this limit. The resulting expressions agree with those 
given in isi[ini[n]. 

For the general case involving a collection of scatterers 
the decay-rate correction functions Fcj and Fc^± show a 
symmetry in the contributions of the electric and mag¬ 
netic multipole amplitudes. This symmetry gets lost for a 
single scatterer. Both correction functions depend on the 
distance of the atom from the center of the spherical 
domain containing the scatterers, and parametrically on 
the radius R of the spherical domain and on the radii a 
of the scatterers themselves (through the scattering am¬ 
plitudes Bf and i?["). These three independent length 
scales that together characterize the configuration (all 
measured in terms of the wavelength of the atomic tran¬ 
sition) occur neatly separated in the summands in (16) 
and (|l7|) . In particular, the dependence on is given by 


^c,ii(^) = +1) E E 


oo oo 


2a3 


i"=c 


xRe 




(16) 


the spherical Hankel functions h^j}\kra) and their deriva¬ 
tives. Since these are proportional to the correc¬ 

tion functions will be oscillating as a function of Va (on 
the scale of the wavelength). For large Va the correction 
functions will decay to 0. 

When the scatterers are taken to be at finite temper¬ 
ature the second moment of the electric field gets an ad¬ 
ditional contribution that is given by the second term in 
(§. The integral over V in that term is taken over all 
positions r' inside the scatterers, so that it is in fact a 
sum over individual contributions for each of the scatter¬ 
ers. As the system of scatterers is dilute, one may use the 
expression (A9) for the Green functions in the integrand. 


F.Ar) = -^ E + 1) E E ^kv.i"{R) 


1=1 


V = 11"=0 


In this way the hot-scatterer contribution to the second 
moment becomes 


xRe 


(-*)*^^ Hi,{kr) 


5{uj — uj')- 


hui^ 


[Ime(a; + iO)] 


E 


dr' 


+(-*)'+! + i?r^r+i'+r') HAkr)] , 


(17) 


7reoc4(e^^‘^ - 1) — j|r'-r,|< 

xGgs(r - ri,r' -ri,u + iO) • Ggg(r - ri,r' -ri,u} + zO). 

( 20 ) 


with the functions F{i{t) = and Hi{t) = 

[d[th[^\t)]/dt]‘^, and with the coefficients 


When the explicit expression (A9) for Ggg is substituted, 
the integration over r' leads to three-dimensional inte¬ 
grals of scala r pro d ucts o f vector spherical wave functions 


cij.m.{R) = {2l'+l){2l"+l)Ii„{R) 


I I' I" 

1 -1 0 


of the form (|All|)-(|A13|). The resulti ng int e grals (a) 


(18) 


may be eliminated with the help of (A15)-(A16). Fur¬ 


thermore, the sum over i may be replaced by an integral 
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over positions r" inside the spherical domain of radius R. 


In this way we get from (20): 
d{u} — u) ) 


oo I 


TTeoc^{eP^ - 1 ) Jv, ^ 

+C^rN[':i(r-r")NS*(r-r") 


1 


2Z + 1 


( 21 ) 


The remaining vector spherical wave functions, which de¬ 
pend on r — r", can be rewritten by using the addition 
theorem of Appendix]^ Finally, we arrive at an expres¬ 
sion for the second moment that is a generalization of 
(15), with the cold-scatterer decay-rate corrections func¬ 


tions Fc^pir) (with p = ||,_L) replaced by new functions 
Fp(r). These contain, apart from the cold-scatterer cor¬ 
rection functions, additional terms that arise when the 
dielectric medium in the scatterers emits thermal radi¬ 
ation. In fact, they have the form Fp{r) = Fc_p(r) -|- 
(e^^“ — 1)“^ Td,p(?’), with the ‘dielectric’ decay-rate cor¬ 
rection functions Fd ||(r) and Fd^j_{r). (The subscript 
d points to the fact that the hot-scatterer contributions 
find their origin in the dielectric medium.) The dielec¬ 
tric decay-rate correction functions have a similar form as 
( 16 ) and ( [T7| ), with the following differences: 1. the func¬ 
tions F[i{t) and Hi{t) are replaced by Fld,i(i) = 
and Hd.iit) = \d\th\^\t)]/dt\'^, respectively; 2. the multi¬ 
pole amplitudes (—Bf (with q = e,m) are replaced 
by — Cf, defined in ( A17|); 3. the symbol Re can be omit¬ 
ted (after the above changes the functions are real). 

The average decay rate (T) for an excited atom in the 
presence of an aggregate of randomly distributed spheri¬ 
cal scatterers at finite temperature gets the form: 


(T) = To + / [A||(r„)ro,|| + Fj.(r„)ro.^] 


( 22 ) 


instead of (19). As we have seen in Section |n| the cold- 
scatterer contribution to the average decay rate follows 
from the first term in which resulted by rewriting 
part of the integral in (1^ by means of the optical the¬ 
orem ([^. Hence, it originates from positions r' both 
in the scatterers and in the surrounding space. In con¬ 
trast, the hot-scatterer contribution has been obtained 
from the second term in ([^, which is an integral over 
the volume V occupied by the dielectric medium within 
the scatterers. To show more clearly the origin of the 
various contributions we may rearrange the total decay- 
rate correction functions Fp by writing Fc^p as the sum 
^d,p + Fr^p of the dielectric decay-rate correction func¬ 
tion Fd^p and a ‘radiative’ decay-rate correction function 
Fi-^p = Fc^p — Fd^p: so that one gets 


^p(^) ~ Fr^p{r) + 




^(3huj _2 


^d.p(r), 


(23) 


with p =11, T. In this form the last term contains all con¬ 
tributions from the dielectric medium within the scat¬ 
terers (both for cold and for hot scatterers), whereas the 


first term represents radiative contributions associated to 
the space outside the scatterers. 

The average absorption rate (Ta) for a ground-state 
atom in the vicinity of a collection of randomly dis¬ 
tributed hot scatterers can likewise be evaluated. As 
absorption and stimulated emission are closely related, 
it is found to be determined by the same correction func¬ 
tions and Fd^j.: 

(ra) = [^’d.||(ra)ro,|| + FdAra)ro,±] . (24) 

The dielectric decay-rate correction functions for hot 
scatterers are quite analogous to their cold-scatterer 
counterparts (16) and In particular, the symmetry 
between electric and magnetic multipole contributions is 
clearly visible in these functions as well. It gets lost in 
the limiting case of a single hot scatterer. As before, the 
independent length scales r^, R and a show up in the 
summand. The dependence on a is contained in the co¬ 
efficients Cf and C[", while R appears in the coefficients 
(18). The distance enters the summand through the 


absolute value of the spherical Hankel functions. As a 
consequence, the oscillating behavior found before is ab¬ 
sent here. In fact, the dielectric correction functions will 
decay smoothly to 0 when the distance increases be¬ 
yond bounds. In the next section we shall determine the 
precise form of the asymptotic behavior for all decay-rate 
correction functions. 


IV. INTEGRAL REPRESENTATIONS FOR THE 
DECAY-RATE CORRECTION FUNCTIONS 

The expressions for the decay-rate correction functions 
as found in the previous section are three-fold sums of 
terms in which the independent variables ra, R and a 
occur in separate factors. To determine their behavior as 
a function of Va, for fixed values of the parameters R and 
a, the sums have to be evaluated numerically. As it turns 
out, a greater numerical efficiency is achieved by starting 
from an alternative representation in which the variables 
Tq and R are intertwined in integral expressions. For cold 
scatterers such a representation can be found by starting 
again from the first term of and using ([Tl])-([T^, as 
before. Upon choosing the atomic position = r = r' 
on the positive z-axis and inserting the representation 
(A8), one may write the zz-component of the integrand 
in (12) as cos^ xGgji {t, u;-|-j0)-|-sin^ xGs,_L(t, w-|-i0), with 


t = — r" I and y the angle between r — r" and the pos¬ 

itive z-axis. The integral over the azimuthal angle p" of 
r" is trivial now. The remaining double integral over the 
spherical variables r" and 9" may be rewritten by intro¬ 
ducing t and u = r" cos 9" as new integration variables, 
so that X is determined by the relation tcosx = r — u. 
The integration over u can be carried out straightfor¬ 
wardly. After insertion of the components Gs,|| and Gs,_l 
of the scattering Green function m. the zz-component 
of (12) leads to an expression for the second moment of 
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the ^-component of the electric field of the same form as Substituting these expressions into (26)-(27l and (31)- 


in (15), with -Fc,||(’’) given as 




/=1 




The functions J®|| ; and J™| j are integrals over t: 


(25) 


rC+P 1 

'^c.im(C) = / dt-{[gi{tX,p) - g2itX,p)]Hi{t) 

JQ-p 


t 




rC+p 

'^c"m(C) = / dtt[gi{tX,p)-g2{tX,p)\Hi{t), 
Jc-p 


(26) 

I, 

(27) 


with p = kR. The auxiliary functions gi are 

P^-{C-t? 


gi{tX,p) = 


g2{t,c,p) = ^ 


2C 

1 


(28) 


2C 


+ 3^ (29) 


while Hi{t) and Hi{t) have been defined below 

Similarly, the a;x-component of the integrand in 


(12) can be written as sin^ y cos^ (/?" Gs_|| (t, w + lO) + 
(cos^ ycos^ (/j" + sin^ (/j")Gs^_L(t, w + iO). Taking the same 
steps as above one arrives at an expression for the sec¬ 
ond moment of the x-component of the electric field as 


in (15). Here, Fc^x(f) is found as 
9 


F.Ar) = - 


16/c^a^ 




1^1 


(30) 


with the integrals: 

'dl±AO = dty\^gi{tx, p) + g2{tx, p)\Hi{t) 

+2l{l+ l)[gi{t,C,,p) - g 2 {t,C,,p)\Hi{t)} , (31) 

/*C+p 

jy\AO= / dtt[gi{t,C„p)+g2{t,C„p)\Hi{t). 

Jc-p 

(32) 

For large R the above results are consistent with those 
found previously m for the correction functions of the 
atomic decay rate in the vicinity of a halfspace. In fact, 
defining = C — p as the (scaled) distance from the 
atom to the surface of the spherical domain, one finds 
the asymptotic forms: 


lim giit,C,p) = t-C 

p—^OO 


1 C'^ 

(33) 


(32) one obtains decay-rate correction functions that are 


in agreement with the results of pT] . 

For hot scatterers similar methods may be used to de¬ 
rive integral representations for F'd,||(^) and Fd,x(^)- To 
that end one starts from (20) and inserts the spherical- 


coordinate representations of the vector spherical wave 
functions in terms of the spherical variables t, y, ip" de¬ 
fined above. Subsequently, one rewrites the triple inte¬ 
gral over r" in terms of the integration variables t, u, p". 
Upon carrying out the integrals over the latter two vari¬ 
ables one gets 






Fd.±{r) = 


(34) 


16 fc 3 ^ 3 E^(^ + l) [CiJl±Akr) 

1^1 

+CrJA±Akr)] , (35) 


with integrals (^°^ P =ll>2- and q = e,m) that 

follow from (26)-(p^ and (31)-(32) by replacing HjJt) 
and Hi{t) with HAil^) and Fld.zUb defined below (21). 


Hence, the expressions for the hot-scatterer integrals are 
quite similar to those for cold scatterers, with squares of 
moduli of spherical Hankel functions (and their deriva¬ 
tives) occurring instead of ordinary squares, as in the 
previous section. 

The integral representations as given above are a suit¬ 
able starting-point to derive the asymptotic behavior of 
the decay-rate correction functions for r tending to oo. 
When r increases, the variable (^ = kr in the Jj-integrals 
(26)-(27) and (31)-(32) gets large. Hence, the integra¬ 


tion variable t in these integrals is large as well, so that 
one may insert the asymptotic form [31] of the spherical 
Hankel function h^At) ~ (—and its derivative 

d[thA\t)]/dt ~ (—i)*e** in the integrands. Instead of t we 
now introduce the new integration variable x, by writing 
t = (^-l-x, so that the integration limits for x become ±p. 
For large one may replace the combinations of gi and 
g 2 in the integrands by the leading terms in their power 
series in 1/^, at fixed values of p and x. For instance, to 
determine the asymptotic form of (26) we write 


l[gi{t,C,p) - g2it,C,p)] - {p^ (36) 


-^g2{t,Q,p)^{p'^-xAl{2A)- 


(37) 


Upon evaluating the (trivial) integral over x we find: 


=2iC 


1 2 3 


'^c,||,z(C) ^ (-1) { ( -7^P + o ) sin(2p) - -pcos(2p) 


-2l{l + l) 


8 


^sin(2p) - ^pcos(2p) 


3 

V 


(38) 
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for large C- In a similar way we get the asymptotic forms 
of the remaining cold-scatterer integrals for large 


JUAO 


o2iC 




--pcos(2p) 




-\p^ + l] sin(2p) 


■ ±,l 


(C)^ 


a2iC 




- sin(2p) — pcos(2p) 


(39) 


.(40) 


and a similar form for J|™(^), which is found to be pro¬ 
portional to instead of . Hence, for increasing 

values of I both integrals get large. However, in the cor¬ 
rection function (25) these integrals are multiplied by the 
multipole amplitudes and summed over all 1. For large I 
the electric-multipole amplitude (as given by (A2)) gets 
the asymptotic form 


Re 


Ime(a;-|-z0) (fca)^*''"^ 

|l-be(w-b*0)|2 22'/2i+2 ’ ^ 


The asymptotic decay-rate correction functions F'c,p(^) 
(with p =|j,T) for cold scatterers, which follow by in¬ 
serting the above asymptotic forms of the Jj-integrals in 
(251 and ( |M| , are damped and oscillating, with a period 
oF(half oifthe wavelength. The damping is algebraic; 
it is proportional to 1/r^ for the longitudinal polariza¬ 
tion, but slower, namely proportional to l/r2, for the 
transverse polarization. For the latter polarization all 
contributions of the electric and magnetic multipoles are 
modulated by the same function of p (or the radius R of 
the spherical domain). Hence, when p is a solution of the 
transcendental equation tan(2p) = 2p, the asymptotic 
form of Fc^±{r) vanishes, which means that the decay is 
faster than l/r2 in that case. Indeed, by evaluating the 
next order in the asymptotic expansion one finds a decay 
proportional to 1/r^ for these special values of R. Hence, 
the scatterers are less effective in modifying the atomic 
decay rate for these particular radii of the spherical do¬ 
main. 

For hot scatterers the asymptotic forms of the J/- 
integrals can be determined in a similar fashion. Since 
only the absolute value of the spherical Hankel functions 
enter the integrands, no oscillations are found. One gets 
the following asymptotic results for large C: 




jm (p, ^ ^ P 

-'d.H.ilU -TF74 


1 

15 


Jl±AO 


(41) 


(42) 


(43) 


The ensuing asymptotic forms of the dielectric decay-rate 
correction functions show a monotonously damped be¬ 
havior proportional to 1/r^ and 1 /r^ for the longitudinal 
and transverse polarizations, respectively. 

For positions near the surface of the spherical domain 
containing the scatterers the decay-rate correction func¬ 
tions turn out to diverge. For cold scatterers and a lon¬ 
gitudinal atomic polarization this follows by considering 
the integrals J®|| j(C) (with q = e,m) for large 1. Since for 
these values the spherical Hankel function gets the form 


h] '{t) ~ i 2*+^/2/^(e^ [3T], the integrals (26) and 


(27) can be evaluated explicitly. One gets: 


■^c.ii,;(C) - 


221 I2l 


- (P ^ -P) 


(44) 


whereas the asymptotic form of the magnetic-multi pole 
amplitude is proportional to l/Z2*+4. Upon inserting (44) 
and (45) into (25) and summing over I one finds that 
the contributions of the electric multipoles diverge for 
C —)■ p—ka, as it becomes a geometric series in ka/{( — p). 
In contrast, the contributions of the magnetic multipoles 
stay finite. As a result we obtain the following asymptotic 
form of the correction function F(,,||(^) for r ^ R + a: 


9 Ime(w -I- iO) R 

^ 16fc3a2 |1 -|- e{uj + z0)|2 [R + a){r — R — a) ’ 

(46) 

so that it diverges linearly near the surface of the spher¬ 
ical domain. Since the scattering spheres may protrude 
from the domain, the actual value of r at which the cor¬ 
rection function diverges is not i?, but i? -|- a. 

The correction function for the transverse atomic po¬ 
larization has a similar divergency, with a prefactor that 
is smaller by a factor of two: Fc^_L{r) ~ These 

asymptotic forms of the decay-rate correction functions 
for an atom near the surface of a spherical domain filled 
with cold scatterers are in agreement with those found 
for an atom near a halfspace with such scatterers [H] , as 
could have been expected. 

Finally, we turn to the asymptotic expressions for the 
decay-rate correction functions pertaining to hot scatter¬ 
ers. Near the surface of the spherical domain containing 
the scatterers these are found to have the same diver¬ 
gency as those for cold scatterers: Fd,||(?’) — -Pc,||(?') and 
Fd,±{r) — Fc^±{r), for r — R + a. The subdominant 
terms in the asymptotic expressions, which are propor¬ 
tional to log(r — i? — a), are found to agree as well. As 
a consequence, the ‘radiative’ decay-rate correction func¬ 
tions Ai.j|(r) and F,-^±{r) stay finite near the surface of 
the domain. 

The linearly divergent behavior of the decay-rate cor¬ 
rection functions and Fdp (with p =||,T) near the 
surface of a spherical domain filled with scatterers is less 
pronounced than that found for the atomic decay rate 
near a single homogeneous sphere [12]. In fact, the latter 
is proportional to l/(r —a)^ for a sphere of radius a. The 
collective effect of an aggregate of scattering spheres is 
mitigated as a consequence of the configurational averag¬ 
ing, which implies a reduced probability for an individual 
scatterer to be near the atom even when is close to 
R Qi. 


[C(C-p)2'+i 




















V. EVALUATION OF THE DECAY-RATE 
CORRECTION FUNCTIONS 


The decay-rate correction functions F^ \\ andT'cT for 


cold scatterers, as given by (25) and (30), depend on 


the distance r via the integrals with p = ||,T and 


q = e,m. Inspection of (26), (27), (31) and (32) shows 


that these integrals can be written as linear combinations 
of a set of basic integrals over the product of two spherical 
Hankel functions with an algebraic prefactor: 

Iiui2,n{z)=J duu~'^hll\u)hll\u), (47) 

for z > 0. As an example, one may write J™| ^ as a linear 
combination of with n = —5, —3, —2, —1,1 and 

z = ( ± p. Clearly, only ‘diagonal’ basic integrals, with 
h = h, show up here. However, the arguments of the in¬ 
tegrals J®|| I and ; in (26) and (31) contain both the 

square of a spherical Hankel function and of the 

derivative (t)]/dt. The latter may be rewritten in 

terms of a spherical Hankel function with a different order 
[ST] , since one has + — 

As a consequence, upon rewriting and J^j_i in 

terms of basic integrals one encounters off-diagonal ba¬ 
sic integrals with I 2 = /i -f 1 as well. It turns out 
that one needs information about the diagonal basic inte¬ 
grals Ippn{z) with n = —5, —3, -2,-1, 0,1,3 and the off- 
diagonal integrals 7/,;+i,n(z) with n = —4, — 2, — 1,0, 2. 
All information about these basic integrals is collected in 
Appendix 

For hot scatterers, the integrals (with p =||,T 

and q = e, m) may be analyzed in a similar way. They 
can be written as linear combinations of a second set 
of basic integrals containing the product of a spherical 
Hankel function and its complex conjugate: 




(48) 


Upon collecting all information about the integrals (47) 


and (48) we can derive explicit expressions for ^ and 


J|p ;, for p = ||,T and q = e,m. Combining these with 
the expressions for the electric and magnetic multipole 
amplitudes (|A2|) we are able to draw curves for the decay- 


rate correction functions. 

For cold scatterers and a transverse atomic polariza¬ 
tion the decay-rate correction function F± = Fc^± is 
shown for a representative choice of R, a and e in Fig. 

As remarked above (23), the correction function Fc^± can 


be split into a dielectric contribution Fd,_L and a radia¬ 
tive contribution F,-^±. Both Fc^± and its constituent 
F). _L show a characteristic damped oscillating behavior, 
which is proportional to for large values of C, in ac¬ 
cordance with the asymptotic form given by (30) with 


(40). In contrast, the dielectric contribution Fd,_L is a 
monotonous function, which falls off like for large 


as has been shown in (43). For atomic positions near the 



FIG. 1: (Color online) Decay-rate correction function F±{() 
(with ^ = kr) for cold scatterers and its constituents Fd,±{() 
(—) and E'_i(C) (-), for a spherical domain with di¬ 

mensionless radius kR — 6 containing cold absorbing spheres 
of dimensionless radius ka = 0.5 and dielectric constant 
e = 3.0 -I- 0.5i. 


surface of the domain containing the scatterers the cor¬ 
rection functions and Fd,j_ diverge as 1/{(—kR—ka), 
in accordance with (46), whereas F^^± stays finite. 

The two contributions to Fc^± are found to be mutually 
counteractive: Fd,j_ is positive, whereas Fr^± is negative. 
Their balance shifts when the dielectric medium becomes 
less absorptive, as is illustrated in Figs. |2|3[ 



FIG. 2: (Color online) Decay-rate correction function F±{() 
for cold scatterers and its constituents Fd,_L(C) (—) and 
Ar,±(C) (-)i for kR = 6, fca = 0.5 and e = 3.0 -I- O.li. 


For hot scatterers the total decay-rate correction func¬ 
tion FI is equal to the sum of Fr _L and — 

l)].Fd,_L’ according to ( |^ . In Fig. an example of the 
total correction function for hot scatterers is represented. 
Comparison with Fig. shows how the balance between 
the dielectric and the radiative contributions gets shifted. 

The decay-rate correction function F’n for the longitu- 
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FIG. 3: (Color online) Decay-rate correction function F±{() FIG. 5: (Color online) Decay-rate correction function f||(C) 
for cold scatterers and its constituents Fd,±(C) (—) and for cold scatterers and its constituents Fd_||(^) (—) and 
^r,±(C) ( -)i for kR = 6, ka = 0.5 and e = 3.0 -I- O.Oli. .F’r.||(C) (-)i for = 6, ka = 0.5 and e = 3.0 -I- 0.5i. 



FIG. 4: (Color online) Total decay-rate correction function 
F±{() for scatterers at finite temperature and its constituents 

[e^"“/(e'^"“-l)]Fd,x(C) (---)andF,,x(C) (-), for kR = 

6, ka — 0.5, £ = 3.0 -I- 0.5i and Phcu = 1. 



FIG. 6: (Color online) Total decay-rate correction function 
1^11 (C) for scatterers at finite temperature and its constituents 

- 1)] Fd.||(0 (- - - ) and F,,,|(C) (-), for kR = 

6, ka — 0.5, £ = 3.0 -I- 0.5i and Phuj = 1. 


dinal polarization direction largely behaves in a similar 
way, as can be seen in Figs. [5]|^ However, in this case 
the oscillatory behavior for large distances shows a faster 
decay, with a damping proportional to in accordance 
with (25) and (38)-(39). The dielectric contribution Fj y 
decays monotonously like as given in (41 )-(42). The 


counterbalancing of the dielectric and radiative contri¬ 
butions occurs here as well, with a shifting balance as 
the absorptive power of the scatterers gets weaker. The 
inversely-linear divergent behavior for small distances, as 
given by (46), is clearly visible. In Fig. it is shown how 
thermal radiation from the scatterers influences the lon¬ 
gitudinal decay-rate correction function: both the total 
correction function and its dielectric part become larger 
owing to stimulated-emission effects, whereas the radia¬ 
tive part remains unchanged. 

The present results for a spherical domain of radius R 


filled with a dilute set of scatterers may be compared to 
those obtained from an effective-medium theory in which 
the domain contains a uniform dielectric with an effec¬ 
tive dielectric constant For scatterers with a size 

much smaller than the wavelength one may use Maxwell 
Garnett theory [53] , with the effective dielectric constant 
j.(eff) = 1 _|_ _ l)/(£ -p 2). For a set of spherical scat¬ 

terers with a radius a that is not small compared to the 
wavelength (as in our model) the effective dielectric con¬ 
stant will depend on ka, as has been discussed in [551135] . 
It may be chosen as = 1 — 3ifBf/{k^a^), with Bf 
the electric-dipole scattering amplitude for a sphere with 
radius a and dielectric constant £, as given in (A2). In 


the limit of small ka the latter expression reduces to the 
Maxwell Garnett form. 

The decay-rate correction functions that follow from 
the effective-medium theory will depend on the electric 
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and magnetic multipole amplitudes (with 

p = e, m) for a uniform sphere with radius R and dielec¬ 
tric constant Since the latter differs from 1 by a 

small amount, we may write the amplitudes R) 

as — l)Bf, with reduced amplitudes Bf depending 

on p = kR: 

~ '^\{i -1-1)^ {(^ +1 + \p'^)[ji(p)\^ + |p^bi-i-i(p)]^ 
-(^ + ^)pji{p)ji+i{p)] , (49) 

= ^'^^^P{\p"\3i{p)f + \p^\ji+i{p)f 
-{l + h)pji{p)ji+i{p)] ■ (50) 



In terms of these reduced amplitudes the effective- 
medium decay-rate correction functions for cold scatter- 
ers become 


W = TO T.W + l)]^Re [Bt{-^)^ 


1=1 


xBf 


kPr'^ 


Hi{kr) 


(51) 


p(eff)/- ^ 9 


KT(^) = 


4k^a^ 


J^l(l + l)Re{B^(-zy 


1^1 


Bf—Jii(kr) + Brffi(kr) 

/c^r^ 


(52) 


The effective-medium dielectric decay-rate correction 
functions T)|'*^^(r’) (with p =||,_L), which come into play 
for hot scatterers, have a similar form. They follow from 
(|^-( 52) by replacing ffi and ffi with the real functions 


Bd,i and ffd,i (defined below (21)) and changing the over¬ 
all sign. It should be noted that (—i^Bi (for q = e,m) 
is real, so that the functions (r) are proportional to 
Re(S(=). 

Having established the expressions for the effective- 
medium decay-rate correction functions, we may com¬ 
pare the ensuing curves with those from our model in 
which the effects of inhomogeneities were taken into ac¬ 
count. As an example, we shall consider the correction 
functions for the choice of parameters considered in Fig. 

From Fig. it is clear that for large distances the 
effective-medium theory gives adequate results for the 
decay-rate correction function F± = Fc,_l for cold scat¬ 
terers, whereas it is off (by about a factor 2) for the 
dielectric part Fd,± and also for the radiative part Fr_j_, 
which play a role for hot scatterers. Thus, the effective- 
medium theory is not able to predict correct values for 
the decay-rate correction functions in the presence of hot 
scatterers. This state of affairs could have been expected 
from the general arguments put forward in [35) . 

In the near zone the predictions of the effective- 
medium theory get even less reliable. As remarked at 
the end of Section |IV[ the small-distance asymptotic be¬ 
havior of the decay-rate correction functions near a do¬ 
main containing a set of discrete scatterers is different 


FIG. 7: (Color online) Comparison of the decay-rate correc¬ 
tion function (—) for a cold effective medium and 

its constituents F^^^bC) (-)> (-) with their 

scattering-medium counterparts F±(Q, Fd,±((), Fr,±((^) (.) 

for cold scatterers, as taken from Fig. 2, for kR = 6, fca = 0.5 
and £ = 3.0 -I- O.li. 


from that for a uniformly-filled domain. Moreover, the 
effective-medium theory misses the contribution of all 
higher-multipole scattering amplitudes, which get impor¬ 
tant near the scatterers. The resulting differences are 
clearly visible in the curves in Fig. 

Further insight in the cause of the discrepancies found 
in Fig.j^is gained by determining the asymptotic behav¬ 
ior of the effective-medium decay-rate correction func¬ 
tions for large kr. Substituting the asymptotic form of 
the spheric al H ankel functions contained in F[i(kr) and 
Hi(kr) in (52), we get an expression for that 

is proportional to ~ Inserting 

(49)-(50), using the rec urren ce relations for ji EH anc. 
employing the sum rule (CIO) we get the asymptotic re¬ 
sult 




sin(2p) - pcos(2p) Re ^ 

(53) 

for large C, = kr. Comparison with (30) and (40) shows 


that we have recovered the contribution of the electric- 
dipole scattering amplitude. Clearly, it dominates the 
long-distance behavior of Fjf*(f^(r), at least for the rather 
small dimensionless radius of the scatterers chosen here 
(ka = 0.5). 

The asymptotic form of Fj®®^(r) may likewise be de¬ 
termined. It is found to be proportional to the sum 
Y^yLiK^ + 1)(—f)^ (A^ + Bjy ), which may be evaluated 
with the help of (C7)-(C9). One gets the asymptotic 
form for large C- 




3 P^ 
2fc3a3 


Re(Ri®) 


(54) 


Comparing with (35) and (43) we see that even the 
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electric-dipole contributions differ: in the effective- 
medium theory the factor Cf (defined in ( |A17 )) is re¬ 
placed by —Re{Bf). Since these differ by a factor 2.06 
in the present case, the discrepancy in Fig. is fully ex¬ 
plained. Incidentally, it should be remarked that choos¬ 
ing a different form for the effective dielectric constant 
(for instance that of the original Maxwell Garnett the¬ 
ory) does not improve matters here. An analysis of the 
decay-rate correction functions for the longitudinal po¬ 
larization (II) and of their asymptotic behavior leads to 
similar conclusions. 


VI. DISCUSSION 

In the previous sections the modification of the decay 
rate of an excited atom in the vicinity of a spherical ag¬ 
gregate of randomly distributed small dielectric spheres 
at arbitrary temperature has been analyzed in detail. 
The changes could be described by a set of decay-rate 
correction functions that depend on three independent 
length scales, namely on the radii of both the scatter- 
ers and of the domain containing them, and on the dis¬ 
tance between the atom and the center of the aggregate. 
Two types of decay-rate correction functions have been 
introduced, each for two different polarizations. The 
decay-rate correction functions and Fc^± describing 
the influence of cold scatterers were shown to exhibit a 
damped oscillatory behavior as a function of the atomic 
distance. In contrast, the dielectric correction functions 
F)j II and Fd._L that come into play for hot scatterers are 
monotonous as a function of the distance. The latter can 
be used as well to split the decay-rate correction functions 
for both cold and hot scatterers in contributions with a 
different origin, viz a dielectric and a radiative part. As 
we have seen, these two parts represent decay-modifying 
effects that may be counteracting, at least for sufficiently 
high absorptive power of the scatterers. 

Two rather different representations for the decay-rate 
correction functions have been found in Sections cni and 


IV The first of these, which has been given in (16) 


(171, follows by using addition theorems for vector spher¬ 


ical wave functions. The resulting expressions are triple 
sums over terms in which the electric and magnetic mul¬ 
tipole contributions show up on an equal footing, while 
the three independent length scales are neatly separated. 
In the second representation, which involves sums of in¬ 
tegrals over finite intervals, two of these length-scales 
(namely, the radius of the domain containing the scat¬ 
terers and the distance of the atom to the center of 
the aggregate) get intertwined in the integrands and the 
boundaries of the integrals. Furthermore, the symme¬ 
try between the electric and magnetic contributions is no 
longer obvious. Nevertheless, this representation is quite 
useful, as it is more suitable for an efficient numerical 
evaluation of the decay-rate correction functions. 


direct consequence of the approximations for the Green 
function that were introduced in 0 and above ( [^ . 
It amounts to the neglect of multiple-scattering events 
that would be important in an aggregate with a high 
filling fraction. When taking into account the effects of 
multiple scattering one expects a nonlinear dependence of 
the decay rate on the density. As an analytical evaluation 
is not feasible in that case, one has to resort to numerical 
simulations in order to treat dense aggregates of spheres. 

A comparison of our results with those obtained 
from an effective-medium theory shows how the inho¬ 
mogeneities brought about by the discrete scatterers af¬ 
fect the atomic decay rate. For scatterers that are far 
from the atom the dielectric correction functions F^.p and 
the radiative correction functions F^^p (for p =T, ||) both 
change considerably when the medium is smoothed, while 
their sums F). p are almost left unchanged. Hence, for 
cold scatterers at large distances from the atom effective- 
medium theory is a useful approximation, whereas it is 
unreliable for scatterers at finite temperature. In con¬ 
trast, for scatterers near the atom the predictions of 
effective-medium theory are inadequate at any tempera¬ 
ture. 

As we have seen, the rate of photon absorption of a 
ground-state atom near a collection of dielectric spheres 
emitting thermal radiation is determined by the same di¬ 
electric decay-rate correction functions F^ ^^ and Fd,_L. If 
the atom is in a mixed state involving its ground state 
and excited state, both atomic decay and atomic ab¬ 
sorption processes may occur, with rates given by (22) 


and (24). In that case a stationary situation may re¬ 
sult with an atomic population ratio Uejug determined 
by the equilibrium condition rig (Ft) = ng{Ta}. Obvi¬ 
ously, for positions far away from the aggregate the atom 
will end in its ground state, since the influence of the 
scatterers gets small. For distances near the surface of 
the domain containing the scatterers all correction func¬ 
tions become large and equal, apart from trivial factors 
2. Hence, the vacuum decay rate Fq in (22) can be ne- 


The changes in the decay rate that follow from (22) are 


proportional to the density Ug of the scatterers. This is a 


glected, so that (F) and (Fa) are proportional in that 
case. As a consequence, the population ratio n^/ng for 
an atom near the surface of the domain reduces to the 
standard Boltzmann factor Glearly, the aggre¬ 

gate of scatterers imposes its temperature on the atom 
in this case. For intermediate atomic positions the pop¬ 
ulation ratio at which emission and absorption processes 
are in equilibrium differs from this simple result, since 
the influence of the scatterers gets smaller with distance. 
The precise form of the population ratio as a function 
of the atomic distance follows directly from the explicit 
expressions for the correction functions that have been 
determined in this paper. 


Appendix A: Green functions 

The Green function in the presence of a dielectric 
sphere, with radius a and complex dielectric constant e, 
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is the sum of the vacuum Green function Gq and a scat¬ 
tering term Gg. For a sphere centered at the origin the 
latter has the form [55H55] 


G.iry,u + ,0) = kf: ■£ 


21 + 1 


Bf + Br M^ir)M^Jr') 


(Al) 


for r and r' both outside the sphere. 

The electric and magnetic multipole amplitudes read 

[miMiiio] 




21 + 1 N[ 

WTi)^^ 


(A2) 


with I > 1 and p = e,m. The numerators and denomi¬ 
nators are given as 


= sfi (<?) ji (?') - ji (q) fi {q') , 

= fM)ji{q') - Mq) /;(?'), 

Df = e (q) ji{q) - (q) fi {q), 

Dr = t\q)jl{q')-h[^\q)Mq'), (A3) 

with fi{q) = {I + l)jiiq) - qji+iiq) and fj;^\q) = {1 + 
l)hi^\q) — qhi^^iq). The spherical Bessel and Hankel 
functions depend on q = ka and q' = \/eq^ with a the 
radius of the sphere. 

The vector spherical wave functions are defined as 


The coefficients with p = e, o are 


=('-11 
,m' ,1" ,m" V 

l V l” \ { l v v 

1 -1 0 




{21 + 1){2V + 1){21" + 1) 


47r 


with symbol s 6f that have been defined below (14). The 
coefhcients (A7) agree with those given in [H], and are 
equivalent with the results obtained in [32] and |43|, as 
follows by employing a couple of identities for 6j-symbols 


The addition theorem for N/„i has the same form as 
( |A6| ), with M and N interchanged. This follows imme¬ 
diately from the identities A and 

Finally, we remark that analogous 
addition theorems hold true for and n|^. In par¬ 


ticular, for r' > r" the addition theorem for M|^(r'-|-r") 
follows from (A6) by replacing both M and N by their 
counterparts and 


For coinciding positions r and r' the scattering term 
in the Green function is diagonal in spherical coor¬ 
dinates: 


Gs(r, r, w + *0) = Ggji (r, uj + * 0 ) 6 ^ 0 ,. 

+Gs,_L(r, w + i0)(eeeg-|-, (A-8) 

with components given in [HIS]. 

If r' is situated inside the sphere and r is outside, the 
full Green function Gq -f Gg = Gqs reads 


M;^„(r) = V A [ripi^^ir)] , (A4) 

N,^„(r) = fc-iy A [V A [rV';.„(r)]], (A5) 


where the scalar spherical wave function stands 

for ji{kr) with ji a spherical Bessel function and 

Yi „i a spherical harmonic depending on the angles {0, ip) 
that determine the direction of the unit vector f = r/r. 
The superscripts (h) in (Al) denote the analogous vector 
spherical wave functions with spherical Hankel functions 


instead of ji. 

The vector spherical wave functions and 

satisfy addition theorems [32| (37] |4T| |l2]. To 
derive them in a concise way |43j one may use an ex¬ 
pansion for I that is a generalization of the standard 
Rayleigh expansion for in terms of spherical har¬ 
monics. The addition theorem for M; „j(r) reads: 


Gos(r,r',w-bi0) =^ 


rn— — l 


2;-bi 




(A9) 


The modified vector spherical wave functions and 

follow from (A4) and (A5) upon replacing k by 
k' = y/ek. Furthermore, the amplitudes are defined 
as 






2/ H- 1 Cp 

I(ITT) 


(AlO) 


oo V oo l" 

M,,„(r' + r") =47r^ E E E 

l' = l m' = -l' l"=0 m"=-l" 


l{l + l) 


/■(®) Mi/ ,(r'^ 

-I- 1) > 


with I > 1 and p = e,m. The coefficients in the numer¬ 
ator are Cg = Cj„ = 1, while the denominators follow 
from (|A3]). 


In the main text we need expressions for integrals over 
scalar products of the modified vector spherical wave 
functions. These can be evaluated by starting from (A4)- 


(A5), employing the standard expressions for the curl of 


(A6) 

a vector in spherical coordinates and evaluating the an- 
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gular integrals. In this way one gets: 


J r<.a 


(All) 


X a + l)/i!Ua) + <\(«)] , (A12) 

[ rfrM|5,(r).N|fA (r)=o, (A13) 

^ r<a 


with the radial integrals(a) = dr r^\ji{k'r)\^. The 
latter are given as: 


ll"\a) = 2a^Re 


k' 


k'^ - {k'*) 


di+i{k'a)ji{k'*a) 


(A14) 


as may be verified by differentiation. This identity is 
useful in proving relations between and Bf. In fact, 
by using the Wronskian for spherical Bessel functions [3T] 
one may derive the following equalities for Z > 1: 


[Ims(cc + ^0)] ll^\a)\AT\^ = (A15) 

[Ime(w + iO)] \{l + |Afp = 


{2[ + lla^ 
l{l + l) k ’■ 


(A16) 


Here we introduced the abbreviations 

Cf = Re , (A17) 


for p = e,m and I > 1. 


Appendix B: Evaluation of integrals 


In the main text we introduced the indefinite integrals 

IiuhA^) = J duu~'^h[l\u)h‘ll\u), (Bl) 

for real z > 0, integer n and non-negative integers li,l 2 - 
We need information about these integrals for the diag¬ 
onal case with li = I 2 = I and for the off-diagonal case 
with Zi = Z 2 — 1 = Z. In the following we shall omit any 
additive constants which may show up when deriving ex¬ 
pressions for these integrals. 

The diagonal integrals with Zi = Z 2 = Z satisfy a recur¬ 
sion relation of the form 


(2Z -|- n -|- 3)//+i^i+i_„(z) — (2Z — u -|- l)/qi^„(z) 


^ —n+1 


hl^\z) 


'b+1 


(^) 


(B2) 


as may be verified by differentiation. When these diago¬ 
nal integrals are known, the non-diagonal integrals follow 
from the relation 


/i.z+i,„(z) = 4(2Z-n)/M.„+i(^)-iz-" ZrW(z) , (B3) 


which likewise may be established by differentiation. 

The recursion relation (B2) can be employed to link 
Ii,i,niz) to a lower-order integral with Iq < 1. 

When n is even, we may choose Zq = 0. For odd n < 
—3 the recursion stops at a positive value of Z. As a 
consequence, one should choose Zq = ^ ® ™ 

that case. With these values of Zq the solution of the 
recursion relation reads for Z > Zg: 


b,/,n 


(z) = 


(Zo - A - 


2^ 


— n+1 


- 2 " 


— n+1 


(^0 H 

i-i 

E 

^—^0 “ 1 “ 1 
(Zq — a + 


1\ 

2/1 —h T 


{2k 


2^-/0 


lQ,lo,n 


A 


{k — in 




(fc + 


' + 2)'- 


fe-l-1 




^-n+l 


(^0 ■ 
1 


2Z -|- n -f 1 


T- + 2 )'-^ 




(B4) 


with Pochhammer symbols (a)„ = a(a-l-l) • • • (a-|-n —I). 
For Z = Zg -I-1 the sum drops out. As we shall see below, 
the sum can be evaluated in closed form for all even n < 
—2 and for all odd n > 1. We shall consider in the 
following the diagonal integrals with even n = 0, — 2 and 
with odd n = 3,1, —1, —3, —5, as these are needed in the 


main text. 

For n = 0 the expression (B4| becomes, upon choosing 
Zg = 0 and evaluating the initial condition in terms of the 
exponential integral as Io,o,o{z) = 2iEi{—2iz) + jz'. 


4i,/,c(z) = — 


2z 

21 + 1 




k^O 


21 + 




2i 


21 + 1 


Ei{—2iz) 


(B5) 


for all Z > 0. The exponential integral can be rewritten 
m in terms of sine and cosine integrals as Ei{—2iz) = 
—Ci(2z) — iSi(2z) -|- iTr/2. 

For n = —2 one gets, again taking Zg = 0 and inserting 
the (trivial) result for Jg,g__ 2 : 


Ii+-2iz) 2(2Z + l)z^ g _ i)(2fc + 3) [4 E)]" 


The sum at the right-hand side can be evaluated with 
the help of (Cl). In that way we arrive at the simpler 
result 


h+-2{z) = ]^z^[h^^\z)f + ^z3[Zi[)J_\(z)]2 
-i(2Z + l)z+|')(z)ZiW(z), (B7) 
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for I > 0. 

Next we consider the diagonal integ rals with odd values 
of n. For n = 3 we choose /q = 0 in (B4). For I > 2 the 


first and third terms at the right-hand side drop out. As 
a result we get: 

= ~ 2{i -\uz^ 

for Z > 2. As before, the sum at the right-hand side can 
be evaluated in explicit form in terms of and 
as is shown in Appendix Upon substituting (C4) for 
p = 1 we arrive at the expression: 




HI - 1)4 


+ 


1 


1 


6(Z-1)/ 

1 


L3(Z-1)4 

2z 


6 (Z-l)(Z + 2 ) 

1 


2 (Z- l)z2_ 

l 2 




3(Z-l)Z(Z-f 2) 3(/-l)z 


(z), (B9) 


for Z > 2. This result is useless for Z = 0,1. These 
special cases can be evaluated in terms of the exponential 


integral, as in (B5|. 


For n = 1 (and Zg = 0) the first term at the right- 
hand side of (B4) drops out for all Z > 1. The special 


case Z = 0, which leads to an exponential integral, has to 
be treated separately. For Z > 1 one gets an expression 
involving a sum that can be evaluated with the help of 
the sum rule (C31. The final result is 


Ii,i,i{z) = 


2Z(Z 


_2l{l + l) 

( 1 ) 


1 

Hi 


[h[^\zW 


'^[hY^iiz)]'^ - , (BIO) 


for Z > 1. 

For n = — 1 we get from (B4) (for Zg = 0) upon insert¬ 
ing the initial condition /o_o,_i(z) = Ei{—2iz): 


Ii,i,-i{z) = -^z^J2 


+ + Ei(-2zz) + , (Bll) 

for all Z > 0. As it turns out, the sum cannot be simplified 
with the help of one of the sum rules in Appendix 
When considering the case n = — 3 in ( |B4[ ) we have 
to choose Zg = 1 and substitute the initial condition for 
Z = 1 that follows by a direct evaluation as Ii^i- 3 {z) = 


2fc-t 1 


k{k + 1) 




This sum may be rewritten with the help of the identity 
(C5) for p = 1. As a result we get for all Z > 0: 


2k+ 1 


+ iz4[Z»«(z)]2 + iz4[Z.W(z)]2 

+ l)Ei{-2iz) 


2iz 


(B12) 


which contains the same sum as (Bll). 

rave to evaluate is 


The last diagonal integral that we 


the case with n = —5. It follows from (B4| by choosing 


Zg = 2 and inserting the initial condition for 1 = 2, which 
has to be calculated by hand as l 2 , 2 ,-b{z) = ^Ei{—2iz) + 
(^iz^ — ^z^—^iz+^)e^^^. One encounters a sum over 
squares of spherical Hankel functions with a coefficient 
containing (k — 2)g in the denominator. Again ( |C5| ) (for 
p = 1, 2) is helpful in reducing this sum, so that we arrive 
at the expression: 


Hi,-5(z) = -^(Z - 1 ) 4 ^" E 


2Zc + l 


+ 

+ 




— (Z - 1)(Z -I- 2) -I- 


^^(^ + 2) + -^^ 


^t^k(k + l)' 


Hihi;_[(z)r 


(Z- l)z^h^^\z)H^{(z) 


( 1 ) 


+ -(Z — 1)4 Ei{—2iz) + — 1)46^*^ 

8 16 


(B13) 


for all Z > 0. The same sum as in (Bll) appears once 
again. 

The off-diagonal integrals Ii,i+i,n{z) with n = 
—4, —2, —1, 0, 2 follow from the above diagonal integrals 
by employing (B3). 


To evaluate the dielectric decay-rate correction func¬ 
tions Fg II and Fg _l we also need information about indef¬ 
inite integrals with spherical Hankel functions and their 
complex conjugates in the integrand: 

Ih,h,niz)=J duu~^hll\u)[h\l\u)]* , (B14) 

for 2: > 0, integer n and non-negative integers hyh- Ex¬ 
plicit expressions for these integrals can be derived in a 
similar way as above. The results are analogous, with 
some subtle differences. For n = 0 and li = I 2 = I for 


instance, one finds on a par with (B5) for all Z > 0: 

i 


Ei{-2iz) + {-\iz + . The ensuing expression for w l^(i )/.'||2 , ^^_|. (1 )/ n |2 

T, . a dm mror QrmaroQ nf dtiormal Har,_ ! 2Z + 1 ^ 4- 1 ^ ' 


Ii^i^-Hz) contains a sum over squares of spherical Han 
kel functions H'^\z), with a coefficient (2Zc + l)/(fc — 1 ) 4 . 


fc =0 


2Z + 1'"' 


(B15) 
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The term with the exponential integral is missing here, 
while the squares of the spherical Hankel functions are 
replaced by the squares of their moduli. We refrain from 
listing here the expressions for the other I[_^ ^ that are 

needed in the evaluation of Td || and _l. 


Appendix C: Sums of squares of spherical Hankel 
functions 

In Appendix sums of squares of spherical Hankel 
functions show up. Some of these may be evaluated as 
simple combinations of and An example is the 
sum occurring in (B6). By induction with respect to I 
one proves the following sum rule for I > 0: 

i 


for Z > 0. The complete hierarchy reads for all p > 0 and 
I > p: 

i 

^(2fc F l){k — p + l)2p[h^k = 

k—p 


E 

U‘=o 


{p — k + l)k{l — p + k + 2)2p-2k ^2k+2 


{p-k+ 




(l)l2 


P 


E 

j 

E 

. fc =0 


,/c=0 
P 


{p — k + 1)k{l — p + k + l)2p-2fc ^2fe+2 
{p-k+ l)k+i 


{p — k + l)fc(Z — p + k + l)2p-2k+l ^2k+l 
{p — k+ l)k+i 




'H 'h+i • 
(C4) 


E 

k=0 


1 


1 


(2fc- l)(2fc 


1 


4(2Z + 1) 




( 1 ) 12 
i+lJ 


_—_t_ \h' 

Sy ^ ~ 4(2Z + 3)^ ' 


( 1)12 


A remainder function independent of I does not occur 
here. 

Apart from these sum rules one may prove several hi¬ 
erarchies of identities relating sums of a similar type as 
above. A first one of these reads as follows 


1 


i 

4^ 


(Cl) 


where we omitted the argument 2 of the spherical Hankel 
functions. It turns out that this sum rule is the first 
in a hierarchy of sum rules with an increasing number 
of factors in the denominator of the coefficient in the 
summand. In fact, one may prove for all p > 0 and all 
I > 0 : 


’ E 

k—p-yi 


+ l r^(l)i2 ^ 

(Al-p)2p+2^ ^ 


2p- 1 

2p 


2k+ 1 

t;ik-P+i)2p^ 


(l)l2 

k J 


1 


-Ah\ 


1 


-Ahi 


E 


2k+ 1 


{k-p - ^)2p+3 


fc=0 2 1 

P 




Id'¥ = 


{p-k + l)fe 


-b 


E 

.k^ 

0 

E 


, fc =0 

p 


fc=0 {p ~ k + ^jk+y^ — p + k + |)2p-2/c+l ^ 

_(p - fc -I- l)fc_I_ 

{p-k + \)k+i{l -pFkF \)2p-2k+i 


1 

2k 


IK 


2p{l—p + 2)2p ^ * 2p(/ — P-b l)2p 

for p > 1 and I > p. The remainder function R 2 ^p{z) 
follows by putting I = p, and discarding the left-hand 
side. In Appendix]^ this identity has been employed for 
p = 1 and p = 2, in order to reduce several sums over 
(i)i 2 squares of spherical Hankel functions with complicated 




( 1 ) 12 
l+l\ 


coefficients to sums with simpler coefficients. 

Yet another hierarchy of identities relates sums with 
coefficients containing an increasing number of factors in 
the numerator: 


(p- Zc -b l)fc 


1 


. fc =0 


(p — fc -b |)fe + i(Z — P -b fc -b |)2p-2fc 




FRi^p{z). 


i+i 

(C2) 


^(2fc-b l)(fc-p-b |)2p-Hl[fcfc T = 


k=0 


The last term, which is independent of I, follows by taking 
Z = 0 on both sides of the identity. For small values of 
p the polynomials in Xjz are of low degree, so that this 
identity yields an efficient way of evaluating the sum, in 
particular for higher Z. 

A second hierarchy of sum rules, with an increasing 
number of factors in the numerator of the summand, can 
be established as well. The first sum rule in this hierarchy 
reads 


Y^(2kF\){hfy = 


2p-b 1 

2(P+1) 


2:2^(2fc-b l)(fc-p-b f)2p-i[4^^]^ 


fc=0 


— P + |)2p-|-l 




— P + |)2p+l 


2(P+1) 

ii-p+y2 p+2^^(i)^(i) 
p -b 1 


2(P+1) 
zhy’h\y^ + R3,p{z ): 


Ahlil 

(C6) 


k=0 


= -z’^lhy'^y - z‘^[h\yy + 2(Z -b l)zhy'’h\y ,(C3) 


for p > 0 and Z > 0. For p = 0 the Pochhammer symbol 
(fc -b |)-i in the sum at the right-hand side should be 
read as l/(fc -b ^), as follows by using the alternative 
form (a)„ = r(a-bn)/r(a) |31j . The remainder function 
R 3 ,p{z) is independent of Z and follows by taking Z = 0. 

Similar sum rules may be established for sums over 
squares of the modulus \hy\z)\ of a spherical Hankel 
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function (with a real argument z), or, even more gen¬ 
erally, for sums over products fi(z)gi{z), with // and gi 
equal to ji, yi, or inde pend ently. For instance, 
one may prove a sum rule like (C2), with all re¬ 
placed by and the product replaced by 

its real part (for real arguments). 

The remainder function is found to be 0 in that case. 
The analogues of the sum rules (lC2|, (C4) and (C5) for 


the modulus of the spherical Hankel functions are useful 
in deriving suitable expressions for the integrals ^ 
in Appendix [B} 

The above sum rules, with squares [ifc(z)]^ instead of 
[h^^\z)Y^ may be employed to derive identities for infi- 
nit e sum s that have been used in Section|Vj The analogue 
of (|C4|) for spherical Bessel functions contains a remain¬ 


der function p!z^P/(|)p. Upon taking the limit I —>• oo, 
the other terms at the right-hand side drop out, so that 
one is left with the equality 


OO 

E 

k—p 


{2k + l){k - p + l)2p[jkiz)]'^ = ■. 

V2^P 


(C7) 


for p > 0. For p = 0 this sum rule is well-known ED, 
while for general p > 0 it agree s w ith an identity in 
[45] . Likewise, the analogue of (C6) for [jk{z)Y con¬ 
tains the remainder function —{—p — ^)2p-i-i[cos(2z) -I- 


{p + sin(2z)/2;]/(2(p -|- 1)). In the limit of infinite I 

the second, third and fourth terms at the right-hand side 
drop out, so that a simple recursion relation connecting 
infinite sums is obtained. It may be combined with the 
initial condition m 


^ (C8) 

k=0 

to derive expressions for all sums of the form 

l){k — p + ^) 2 p+i[jk{z)]'^ with p > 0. In particular, one 

gets for p = 0: 

^(2fc -h l)2[jfc(z)]2 = zSi{2z) + ^ cos(22:) , 

(C9) 

which is consistent with an identity involving the gener¬ 
alized hypergeometric function iF 2 (i; |, in [45] . 

Finally, an infinite sum rule with alternating signs: 

OO 

'^{-l)’'{2k + l){k-p+l)2p[jk{z)]'^ = {-lfp\zPjp{2z) 

k—p 

(CIO) 

has been employed in Section |V] for p = 0,1. For p = 0 
it is well-known m, whereas for general p > 0 it can be 
found in [45] . 
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